home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / srcuc.zip / IMAGE.C < prev    next >
C/C++ Source or Header  |  1991-12-20  |  38KB  |  1,234 lines
  1. /* -*-C-*-
  2.  
  3. $Header: /scheme/users/cph/src/microcode/RCS/image.c,v 9.32 1991/12/20 22:49:23 cph Exp $
  4.  
  5. Copyright (c) 1987-91 Massachusetts Institute of Technology
  6.  
  7. This material was developed by the Scheme project at the Massachusetts
  8. Institute of Technology, Department of Electrical Engineering and
  9. Computer Science.  Permission to copy this software, to redistribute
  10. it, and to use it for any purpose is granted, subject to the following
  11. restrictions and understandings.
  12.  
  13. 1. Any copy made of this software must include this copyright notice
  14. in full.
  15.  
  16. 2. Users of this software agree to make their best efforts (a) to
  17. return to the MIT Scheme project any improvements or extensions that
  18. they make, so that these may be included in future releases; and (b)
  19. to inform MIT of noteworthy uses of this software.
  20.  
  21. 3. All materials developed as a consequence of the use of this
  22. software shall duly acknowledge such use, in accordance with the usual
  23. standards of acknowledging credit in academic research.
  24.  
  25. 4. MIT has made no warrantee or representation that the operation of
  26. this software will be error-free, and MIT is under no obligation to
  27. provide any services, by way of maintenance, update, or otherwise.
  28.  
  29. 5. In conjunction with products arising from the use of this material,
  30. there shall be no use of the name of the Massachusetts Institute of
  31. Technology nor of any adaptation thereof in any advertising,
  32. promotional, or sales literature without prior written consent from
  33. MIT in each case. */
  34.  
  35. #include "scheme.h"
  36. #include "prims.h"
  37. #include "array.h"
  38. #include <math.h>
  39.  
  40. void
  41. arg_image (arg_number, nrows, ncols, array)
  42.      int arg_number;
  43.      long * nrows;
  44.      long * ncols;
  45.      REAL ** array;
  46. {
  47.   fast SCHEME_OBJECT argument = (ARG_REF (arg_number));
  48.   fast SCHEME_OBJECT rest;
  49.   fast SCHEME_OBJECT first;
  50.   fast SCHEME_OBJECT second;
  51.   fast SCHEME_OBJECT third;
  52.   if (! (PAIR_P (argument))) goto loser;
  53.   first = (PAIR_CAR (argument));
  54.   if (! (UNSIGNED_FIXNUM_P (first))) goto loser;
  55.   rest = (PAIR_CDR (argument));
  56.   if (! (PAIR_P (rest))) goto loser;
  57.   second = (PAIR_CAR (rest));
  58.   if (! (UNSIGNED_FIXNUM_P (second))) goto loser;
  59.   rest = (PAIR_CDR (rest));
  60.   if (! (PAIR_P (rest))) goto loser;
  61.   third = (PAIR_CAR (rest));
  62.   if (! (ARRAY_P (third))) goto loser;
  63.   if ((PAIR_CDR (rest)) != EMPTY_LIST) goto loser;
  64.   (*nrows) = (UNSIGNED_FIXNUM_TO_LONG (first));
  65.   (*ncols) = (UNSIGNED_FIXNUM_TO_LONG (second));
  66.   (*array) = (ARRAY_CONTENTS (third));
  67.   return;
  68.  loser:
  69.   error_bad_range_arg (arg_number);
  70.   /* NOTREACHED */
  71. }
  72.  
  73. #define MAKE_IMAGE(nrows, ncols, array)                    \
  74.   (cons ((LONG_TO_UNSIGNED_FIXNUM (nrows)),                \
  75.      (cons ((LONG_TO_UNSIGNED_FIXNUM (ncols)),            \
  76.         (cons ((array), EMPTY_LIST))))))
  77.  
  78. static int
  79. read_byte (fp)
  80.      fast FILE * fp;
  81. {
  82.   int result = (getc (fp));
  83.   if (ferror (fp))
  84.     error_external_return ();
  85.   return (result);
  86. }
  87.  
  88. static int
  89. read_word (fp)
  90.      fast FILE * fp;
  91. {
  92.   int result = (getw (fp));
  93.   if (ferror (fp))
  94.     error_external_return ();
  95.   return (result);
  96. }
  97.  
  98. static void
  99. write_word (fp, datum)
  100.      fast FILE * fp;
  101.      int datum;
  102. {
  103.   if ((putw (datum, fp)) != 0)
  104.     error_external_return ();
  105.   return;
  106. }
  107.  
  108. static int
  109. read_2bint (fp)
  110.      fast FILE * fp;
  111. {
  112.   int msd = (getc (fp));
  113.   if (ferror (fp))
  114.     error_external_return ();
  115.   {
  116.     int lsd = (getc (fp));
  117.     if (ferror (fp))
  118.       error_external_return ();
  119.     {
  120.       int result = ((msd << 8) | lsd);
  121.       return (((result & (1 << 15)) == 0) ? result : ((-1 << 16) | result));
  122.     }
  123.   }
  124. }
  125.  
  126. static void
  127. write_2bint (fp, datum)
  128.      fast FILE * fp;
  129.      int datum;
  130. {
  131.   if (((putc (((datum >> 8) & 0xFF), fp)) == EOF) ||
  132.       ((putc ((datum & 0xFF), fp)) == EOF))
  133.     error_external_return ();
  134.   return;
  135. }
  136.  
  137.  
  138. DEFINE_PRIMITIVE ("IMAGE-READ-ASCII", Prim_read_image_ascii, 1, 1, 0)
  139. {
  140.   fast FILE * fp;
  141.   long nrows, ncols;
  142.   PRIMITIVE_HEADER (1);
  143.   CHECK_ARG (1, STRING_P);
  144.   if ( (fp = fopen((STRING_ARG (1)), "r")) == NULL) 
  145.     error_bad_range_arg (1);
  146.   fscanf (fp, "%d %d \n", (&nrows), (&ncols));
  147.   if ((ferror (fp)) || ((ncols > 512) || (nrows > 512)))
  148.     { printf("read-image-ascii-file: problem with rows,cols \n");
  149.       error_bad_range_arg (1); }
  150.   {
  151.     fast long length = (nrows * ncols);
  152.     SCHEME_OBJECT array = (allocate_array (length));
  153.     fast REAL * scan = (ARRAY_CONTENTS (array));
  154.     while ((length--) > 0)
  155.       { long number;
  156.     fscanf (fp, "%d", (&number));
  157.     if (ferror (fp))
  158.       error_external_return ();
  159.     (*scan++) = ((REAL) number);
  160.       }
  161.     if ((fclose (fp)) != 0) error_external_return ();
  162.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  163.   }
  164. }
  165.  
  166. DEFINE_PRIMITIVE ("IMAGE-READ-2BINT", Prim_read_image_2bint, 1, 1, 0)
  167. {
  168.   FILE *fp, *fopen();
  169.   PRIMITIVE_HEADER (1);
  170.   CHECK_ARG (1, STRING_P);
  171.   if ( ( fp = (fopen((STRING_ARG (1)), "r")) ) == NULL) 
  172.     error_bad_range_arg (1);
  173.   {
  174.     int nrows = (read_word (fp));
  175.     int ncols = (read_word (fp));
  176.     fast long length = (nrows * ncols);
  177.     SCHEME_OBJECT array = (allocate_array (length));
  178.     fast REAL * scan = (ARRAY_CONTENTS (array));
  179.     while ((length--) > 0)
  180.       (*scan++) = ((REAL) (read_2bint (fp)));
  181.     if ((fclose (fp)) != 0)
  182.       error_external_return ();
  183.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  184.   }
  185. }
  186.  
  187. DEFINE_PRIMITIVE ("IMAGE-READ-CTSCAN", Prim_read_image_ctscan, 1, 1, 0)
  188. {
  189.   fast FILE * fp;
  190.   PRIMITIVE_HEADER (1);
  191.   CHECK_ARG (1, STRING_P);
  192.   fp = (fopen((STRING_ARG (1)), "r"));
  193.   if (fp == ((FILE *) 0))
  194.     error_bad_range_arg (1);
  195.   Primitive_GC_If_Needed (BYTES_TO_WORDS (512 * (sizeof (int))));
  196.   {
  197.     int nrows = 512;
  198.     int ncols = 512;
  199.     SCHEME_OBJECT array = (allocate_array (nrows * ncols));
  200.     REAL * Array = (ARRAY_CONTENTS (array));
  201.     fast int * Widths = ((int *) Free);
  202.     fast int i;
  203.     /* Discard header */
  204.     for (i = 0; (i < 2048); i += 1)
  205.       (void) read_byte (fp);
  206.     for (i = 0; (i < 512); i += 1)
  207.       (Widths [i]) = (read_2bint (fp));
  208.     for (i = 0; (i < (nrows * ncols)); i += 1)
  209.       (Array [i]) = 0;
  210.     for (i = 0; (i < 512); i += 1)
  211.       {
  212.     fast int array_index = ((i * 512) + (256 - (Widths [i])));
  213.     fast int m;
  214.     for (m = 0; (m < (2 * (Widths [i]))); m += 1)
  215.       (Array [array_index + m]) = ((REAL) (read_2bint (fp)));
  216.       }
  217.     /* CTSCAN images are upside down */
  218. #if (REAL_IS_DEFINED_DOUBLE != 0)
  219.     ALIGN_FLOAT (Free);
  220.     Free += 1;
  221. #endif
  222.     Primitive_GC_If_Needed (512 * REAL_SIZE);
  223.     Image_Mirror_Upside_Down (Array, nrows, ncols, ((REAL *) Free));
  224.     if ((fclose (fp)) != 0)
  225.       error_external_return ();
  226.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  227.   }
  228. }
  229.  
  230.  
  231. DEFINE_PRIMITIVE ("IMAGE-READ-CBIN", Prim_read_image_cbin, 1, 1, 0)
  232. {
  233.   fast FILE * fp;
  234.   PRIMITIVE_HEADER (1);
  235.   CHECK_ARG (1, STRING_P);
  236.   fp = (fopen ((STRING_ARG (1)), "r"));
  237.   if (fp == ((FILE *) 0))
  238.     error_bad_range_arg (1);
  239.   {
  240.     int nrows = (read_word (fp));
  241.     int ncols = (read_word (fp));
  242.     long length = (nrows * ncols);
  243.     SCHEME_OBJECT array = (allocate_array (length));
  244.     fast REAL * scan = (ARRAY_CONTENTS (array));
  245.     while ((length--) > 0)
  246.       (*scan++) = (read_word (fp));
  247.     if ((fclose (fp)) != 0)
  248.       error_external_return ();
  249.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  250.   }
  251. }
  252.  
  253.  
  254. Image_Mirror_Upside_Down (Array,nrows,ncols,Temp_Row)
  255.      REAL * Array;
  256.      long nrows;
  257.      long ncols;
  258.      REAL * Temp_Row;
  259. { int i;
  260.   REAL *M_row, *N_row;
  261.   for (i=0;i<(nrows/2);i++) {
  262.     M_row = Array + (i * ncols);
  263.     N_row = Array + (((nrows-1)-i) * ncols);
  264.     C_Array_Copy(N_row,    Temp_Row, ncols);
  265.     C_Array_Copy(M_row,    N_row,    ncols);
  266.     C_Array_Copy(Temp_Row, M_row,    ncols);
  267.   }
  268. }
  269.  
  270. /* The following does not work, to be fixed. */
  271. DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0)
  272. {
  273.   long Length;
  274.   long i,j;
  275.   long allocated_cells;
  276.   long nrows, ncols;
  277.   SCHEME_OBJECT array;
  278.   double *Array, *From_Here;
  279.   fast double temp_value_cell;
  280.   float *To_Here;
  281.   PRIMITIVE_HEADER (1);
  282.   arg_image (1, (&nrows), (&ncols), (&array));
  283.   Array = ((double *) (ARRAY_CONTENTS (array)));
  284.   From_Here = Array;
  285.   To_Here = ((float *) (Array));
  286.   Length = nrows * ncols;
  287.  
  288.   for (i=0;i<Length;i++) {
  289.     temp_value_cell = *From_Here;
  290.     From_Here++;
  291.     *To_Here = ((float) temp_value_cell);
  292.     To_Here++;
  293.   }
  294.   /* and now SIDE-EFFECT the ARRAY_HEADER */
  295.   SET_VECTOR_LENGTH (array, ((Length * FLOAT_SIZE) + 1));
  296.   PRIMITIVE_RETURN (UNSPECIFIC);
  297. }
  298.  
  299.  
  300. DEFINE_PRIMITIVE ("SUBIMAGE-COPY!", Prim_subimage_copy, 12, 12, 0)
  301. {
  302.   long r1, c1, r2, c2, at1r, at1c, at2r, at2c, mr, mc;
  303.   REAL * x;
  304.   REAL * y;
  305.   void subimage_copy ();
  306.   PRIMITIVE_HEADER (12);
  307.   CHECK_ARG (3, ARRAY_P);    /* image array 1 = source array */
  308.   CHECK_ARG (6, ARRAY_P);    /* image array 2 = destination array */
  309.   r1 = (arg_nonnegative_integer (1));
  310.   c1 = (arg_nonnegative_integer (2));
  311.   if ((r1 * c1) != (ARRAY_LENGTH (ARG_REF (3))))
  312.     error_bad_range_arg (3);
  313.   x = (ARRAY_CONTENTS (ARG_REF (3)));
  314.   r2 = arg_nonnegative_integer (4);
  315.   c2 = arg_nonnegative_integer (5);
  316.   if ((r2 * c2) != (ARRAY_LENGTH (ARG_REF (6))))
  317.     error_bad_range_arg (6);
  318.   y = (ARRAY_CONTENTS (ARG_REF (6)));
  319.   at1r = (arg_nonnegative_integer (7));
  320.   at1c = (arg_nonnegative_integer (8));
  321.   at2r = (arg_nonnegative_integer (9));
  322.   at2c = (arg_nonnegative_integer (10));
  323.   mr = (arg_nonnegative_integer (11));
  324.   mc = (arg_nonnegative_integer (12));
  325.   if (((at1r + mr) > r1) || ((at1c + mc) > c1))
  326.     error_bad_range_arg (7);
  327.   if (((at2r + mr) > r2) || ((at2c + mc) > c2))
  328.     error_bad_range_arg (9);
  329.   subimage_copy (x, y, r1, c1, r2, c2, at1r, at1c, at2r, at2c, mr, mc);
  330.   PRIMITIVE_RETURN (UNSPECIFIC);
  331. }
  332.  
  333. void
  334. subimage_copy (x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc)
  335.      REAL *x,*y;
  336.      long r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc;
  337. { long i,j;
  338.   REAL *xrow,*yrow;
  339.  
  340.   xrow = x + at1r*c1   + at1c;
  341.   yrow = y + at2r*c2   + at2c;    /*  A(i,j)--->Array[i*ncols+j]  */
  342.  
  343.   for (i=0; i<mr; i++) {
  344.     for (j=0; j<mc; j++)    yrow[j] = xrow[j];
  345.     xrow = xrow + c1;
  346.     yrow = yrow + c2;
  347.   }
  348. }
  349.  
  350. /* image-operation-2
  351.    groups together procedures     that use 2 image-arrays
  352.    (usually side-effecting the 2nd image, but not necessarily) */
  353.  
  354. DEFINE_PRIMITIVE ("IMAGE-OPERATION-2!", Prim_image_operation_2, 5,5, 0)
  355. { long rows, cols, nn, opcode;
  356.   REAL *x,*y;
  357.   void image_laplacian();
  358.   PRIMITIVE_HEADER (5);
  359.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  360.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  361.   CHECK_ARG (3, FIXNUM_P);    /* cols */
  362.   CHECK_ARG (4, ARRAY_P);    /* image array 1 */
  363.   CHECK_ARG (5, ARRAY_P);    /* image array 2 */
  364.  
  365.   x = ARRAY_CONTENTS(ARG_REF(4));
  366.   y = ARRAY_CONTENTS(ARG_REF(5));
  367.   rows = arg_nonnegative_integer(2);
  368.   cols = arg_nonnegative_integer(3);
  369.   nn = rows*cols;
  370.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  371.   if (nn != ARRAY_LENGTH(ARG_REF(5)))   error_bad_range_arg(5);
  372.  
  373.   opcode = arg_nonnegative_integer(1);
  374.  
  375.   if (opcode==1)
  376.     image_laplacian(x,y,rows,cols); /* result in y */
  377.   else if (opcode==2)
  378.     error_bad_range_arg(1);    /* illegal opcode */
  379.   else
  380.     error_bad_range_arg(1);    /* illegal opcode */
  381.  
  382.   PRIMITIVE_RETURN (UNSPECIFIC);
  383. }
  384.  
  385. /* Laplacian form [4,-1,-1,-1,-1]/4
  386.    A(i,j) --> Array[i*ncols + j]
  387.    With no knowledge outside boundary,
  388.    assume laplace(edge-point)=0.0 (no wrap-around, no artificial bndry) */
  389. void
  390. image_laplacian (x,y, nrows,ncols)
  391.      REAL *x, *y;
  392.      long nrows, ncols;
  393. { long i,j, nrows1, ncols1;
  394.   nrows1=nrows-1; ncols1=ncols-1;
  395.   /* no need todo anything for 1-point image */
  396.   if ((nrows<2)||(ncols<2))
  397.     return;
  398.   /* */
  399.   i=0;j=0;           y[i*ncols+j] = 0.0; /* NE corner */
  400.   i=0;j=ncols1;      y[i*ncols+j] = 0.0; /* NW corner */
  401.   i=nrows1;j=0;      y[i*ncols+j] = 0.0; /* SE corner */
  402.   i=nrows1;j=ncols1; y[i*ncols+j] = 0.0; /* SW corner */
  403.   i=0; for (j=1;j<ncols1;j++)       y[i*ncols+j] = 0.0;    /* NORTH row */
  404.   i=nrows1; for (j=1;j<ncols1;j++)  y[i*ncols+j] = 0.0;    /* SOUTH row */
  405.   j=0; for (i=1;i<nrows1;i++)       y[i*ncols+j] = 0.0;    /* EAST column */
  406.   j=ncols1; for (i=1;i<nrows1;i++)  y[i*ncols+j] = 0.0;    /* WEST column */
  407.   /* */
  408.   for (i=1;i<nrows1;i++)
  409.     for (j=1;j<ncols1;j++)
  410.       y[i*ncols+j] =
  411.     x[i*ncols+j] -
  412.       ((x[i*ncols+(j-1)] +
  413.         x[i*ncols+(j+1)] +
  414.         x[(i-1)*ncols+j] +
  415.         x[(i+1)*ncols+j])
  416.        / 4);
  417.   return;
  418. }
  419.  
  420. DEFINE_PRIMITIVE ("IMAGE-DOUBLE-BY-INTERPOLATION", Prim_image_double_by_interpolation, 1, 1, 0)
  421. { long nrows, ncols, Length;
  422.   SCHEME_OBJECT Parray;
  423.   REAL *Array, *To_Here;
  424.   SCHEME_OBJECT Result, array;
  425.   PRIMITIVE_HEADER (1);
  426.   arg_image (1, (&nrows), (&ncols), (&Parray));
  427.   Array = (ARRAY_CONTENTS (Parray));
  428.   array = (allocate_array (4 * nrows * ncols));
  429.   Result = (MAKE_IMAGE (nrows, ncols, array));
  430.  
  431.   Array = ARRAY_CONTENTS(Parray);
  432.   C_image_double_by_interpolation
  433.     (Array, (ARRAY_CONTENTS(array)), nrows, ncols);
  434.   PRIMITIVE_RETURN(Result);
  435. }
  436.  
  437. /* double image by linear interpolation.
  438.    ---new_array must be 4 times as long ---
  439.    A(i,j) --> Array[i*ncols + j]
  440.    magnification in a south-east direction
  441.    (i.e. replication of pixels in South-East corner) */
  442. C_image_double_by_interpolation (array, new_array, nrows, ncols)
  443.      REAL *array, *new_array;
  444.      long nrows, ncols;
  445. { long i,j, nrows1, ncols1, nrows2, ncols2;
  446.   nrows1=nrows-1; ncols1=ncols-1;
  447.   nrows2=2*nrows; ncols2=2*ncols;
  448.   /* no need todo anything for 1-point image */
  449.   if ((nrows<2)||(ncols<2)) return(1);
  450.   i=nrows1; for (j=0;j<ncols1;j++)     /* SOUTH row */
  451.   { new_array[(2*i)*ncols2+(2*j)]      = array[i*ncols+j];
  452.     new_array[(2*i+1)*ncols2+(2*j)]    = array[i*ncols+j];
  453.     new_array[(2*i)*ncols2+(2*j)+1] =
  454.       ((array[i*ncols+j]+array[i*ncols+j+1]) / 2);
  455.     new_array[(2*i+1)*ncols2+(2*j)+1]  = new_array[(2*i)*ncols2+(2*j)+1];
  456.   }
  457.   j=ncols1; for (i=0;i<nrows1;i++)      /* WEST column */
  458.   { new_array[(2*i)*ncols2+(2*j)]      = array[i*ncols+j];
  459.     new_array[(2*i)*ncols2+(2*j)+1]    = array[i*ncols+j];
  460.     new_array[(2*i+1)*ncols2+(2*j)] =
  461.       ((array[i*ncols+j]+array[(i+1)*ncols+j]) / 2);
  462.     new_array[(2*i+1)*ncols2+(2*j)+1]  = new_array[(2*i+1)*ncols2+(2*j)];
  463.   }
  464.   i=nrows1;j=ncols1; {                  /* SW corner */
  465.     new_array[(2*i)*ncols2+(2*j)]     =  array[i*ncols+j];
  466.     new_array[(2*i)*ncols2+(2*j)+1]   =  array[i*ncols+j];
  467.     new_array[(2*i+1)*ncols2+(2*j)]   =  array[i*ncols+j];
  468.     new_array[(2*i+1)*ncols2+(2*j)+1] =  array[i*ncols+j];
  469.   }
  470.   /* */
  471.   for (i=0;i<nrows1;i++)
  472.     for (j=0;j<ncols1;j++) {                        /* interior of image */
  473.       new_array[(2*i)*ncols2+(2*j)] =  array[i*ncols+j];
  474.       new_array[(2*i)*ncols2+(2*j)+1] =
  475.     ((array[i*ncols+j]+array[i*ncols+j+1]) / 2);
  476.       new_array[(2*i+1)*ncols2+(2*j)] =
  477.     ((array[i*ncols+j]+array[(i+1)*ncols+j]) / 2);
  478.       new_array[(2*i+1)*ncols2+(2*j)+1] =
  479.     ((array[i*ncols+j] +
  480.       array[i*ncols+j+1] +
  481.       array[(i+1)*ncols+j] +
  482.       array[(i+1)*ncols+j+1])
  483.      / 4);
  484.     }
  485. }
  486.  
  487. DEFINE_PRIMITIVE ("IMAGE-MAKE-RING", Prim_image_make_ring, 4, 4, 0)
  488. {
  489.   PRIMITIVE_HEADER (4);
  490.   {
  491.     long nrows = (arg_nonnegative_integer (1));
  492.     long ncols = (arg_nonnegative_integer (2));
  493.     long Length = (nrows * ncols);
  494.     long high_cycle =
  495.       (arg_index_integer (4, ((min ((nrows / 2), (ncols / 2))) + 1)));
  496.     long low_cycle = (arg_index_integer (3, (high_cycle + 1)));
  497.     SCHEME_OBJECT Ring_Array_Result = (allocate_array (Length));
  498.     SCHEME_OBJECT Result = (MAKE_IMAGE (nrows, ncols, Ring_Array_Result));
  499.     REAL * Ring_Array = (ARRAY_CONTENTS (Ring_Array_Result));
  500.     C_Image_Make_Ring (Ring_Array, nrows, ncols, low_cycle, high_cycle);
  501.     PRIMITIVE_RETURN (Result);
  502.   }
  503. }
  504.  
  505. C_Image_Make_Ring (Ring_Array, nrows, ncols, low_cycle, high_cycle)
  506.      REAL *Ring_Array;
  507.      long nrows, ncols, low_cycle, high_cycle;
  508. { long Square_LC=low_cycle*low_cycle, Square_HC=high_cycle*high_cycle;
  509.   long i, j, m, n, radial_cycle;
  510.   long nrows2=nrows/2, ncols2=ncols/2;
  511.   for (i=0; i<nrows; i++) {
  512.     for (j=0; j<ncols; j++) {
  513.       m = ((i<nrows2) ? i : (nrows-i));
  514.       n = ((j<ncols2) ? j : (ncols-j));
  515.       radial_cycle = (m*m)+(n*n);
  516.       if ( (radial_cycle<Square_LC) || (radial_cycle>Square_HC))
  517.     Ring_Array[i*ncols+j] = 0;
  518.       else Ring_Array[i*ncols+j] = 1;
  519.     }}
  520. }
  521.  
  522. /* Periodic-shift without side-effects for code simplicity. */
  523.  
  524. DEFINE_PRIMITIVE ("IMAGE-PERIODIC-SHIFT", Prim_image_periodic_shift, 3, 3, 0)
  525. {
  526.   long nrows;
  527.   long ncols;
  528.   REAL * Array;
  529.   PRIMITIVE_HEADER (3);
  530.   {
  531.     SCHEME_OBJECT Parray;
  532.     arg_image (1, (&nrows), (&ncols), (&Parray));
  533.     Array = (ARRAY_CONTENTS (Parray));
  534.   }
  535.   {
  536.     long ver_shift = ((arg_integer (2)) % nrows);
  537.     long hor_shift = ((arg_integer (3)) % ncols);
  538.     SCHEME_OBJECT array = (allocate_array (nrows * ncols));
  539.     SCHEME_OBJECT Result = (MAKE_IMAGE (nrows, ncols, array));
  540.     C_Image_Periodic_Shift
  541.       (Array,
  542.        (ARRAY_CONTENTS (array)),
  543.        nrows,
  544.        ncols,
  545.        ver_shift,
  546.        hor_shift);
  547.     PRIMITIVE_RETURN (Result);
  548.   }
  549. }
  550.  
  551. /* ASSUMES ((hor_shift < nrows) && (ver_shift < ncols)) */
  552.  
  553. C_Image_Periodic_Shift(Array, New_Array, nrows, ncols, ver_shift, hor_shift)
  554.      REAL *Array, *New_Array; long nrows, ncols, hor_shift, ver_shift;
  555. { long i, j, ver_index, hor_index;
  556.   REAL *To_Here;
  557.   To_Here = New_Array;
  558.   for (i=0;i<nrows;i++) {
  559.     for (j=0;j<ncols;j++) {
  560.       ver_index = (i+ver_shift) % nrows;
  561.       if (ver_index<0) ver_index = nrows+ver_index; /* wrapping around */
  562.       hor_index = (j+hor_shift) % ncols;
  563.       if (hor_index<0) hor_index = ncols+hor_index;
  564.       *To_Here++ = Array[ver_index*ncols + hor_index];
  565.     }}
  566. }
  567.  
  568. /* Rotations and stuff */
  569.  
  570. DEFINE_PRIMITIVE ("IMAGE-TRANSPOSE!", Prim_image_transpose, 4,4, 0)
  571. { long rows, cols;
  572.   REAL *x, *y;
  573.   PRIMITIVE_HEADER (4);
  574.   CHECK_ARG (1, FIXNUM_P);    /* rows */
  575.   CHECK_ARG (2, FIXNUM_P);    /* cols */
  576.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  577.   CHECK_ARG (4, ARRAY_P);    /* image array 2,  empty for rows=cols */
  578.   
  579.   rows = arg_nonnegative_integer(1);
  580.   cols = arg_nonnegative_integer(2);
  581.   x = (ARRAY_CONTENTS (ARG_REF(3)));
  582.   y = (ARRAY_CONTENTS (ARG_REF(4)));
  583.   
  584.   if (rows==cols)        /* square image ==> ignore argument 4  */
  585.     Image_Fast_Transpose (x, rows);
  586.   else
  587.     Image_Transpose (x, y, rows, cols);
  588.   
  589.   PRIMITIVE_RETURN (UNSPECIFIC);
  590. }
  591.  
  592. DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CLW!", Prim_image_rotate_90clw, 1, 1, 0)
  593. {
  594.   long nrows;
  595.   long ncols;
  596.   REAL * Array;
  597.   long Length;
  598.   PRIMITIVE_HEADER (1);
  599.   {
  600.     SCHEME_OBJECT Parray;
  601.     arg_image (1, (&nrows), (&ncols), (&Parray));
  602.     Array = (ARRAY_CONTENTS (Parray));
  603.   }
  604.   Length = (nrows * ncols);
  605. #if (REAL_IS_DEFINED_DOUBLE != 0)
  606.     ALIGN_FLOAT (Free);
  607.     Free += 1;
  608. #endif
  609.   Primitive_GC_If_Needed (Length * REAL_SIZE);
  610.   {
  611.     REAL * Temp_Array = ((REAL *) Free);
  612.     Image_Rotate_90clw (Array, Temp_Array, nrows, ncols);
  613.     C_Array_Copy (Temp_Array, Array, Length);
  614.   }
  615.   {
  616.     SCHEME_OBJECT argument = (ARG_REF (1));
  617.     SET_PAIR_CAR (argument, (LONG_TO_UNSIGNED_FIXNUM (ncols)));
  618.     SET_PAIR_CAR ((PAIR_CDR (argument)), (LONG_TO_UNSIGNED_FIXNUM (nrows)));
  619.     PRIMITIVE_RETURN (argument);
  620.   }
  621. }
  622.  
  623. DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CCLW!", Prim_image_rotate_90cclw, 1, 1, 0)
  624. {
  625.   long nrows;
  626.   long ncols;
  627.   REAL * Array;
  628.   long Length;
  629.   PRIMITIVE_HEADER (1);
  630.   {
  631.     SCHEME_OBJECT Parray;
  632.     arg_image (1, (&nrows), (&ncols), (&Parray));
  633.     Array = (ARRAY_CONTENTS (Parray));
  634.   }
  635.   Length = (nrows * ncols);
  636. #if (REAL_IS_DEFINED_DOUBLE != 0)
  637.     ALIGN_FLOAT (Free);
  638.     Free += 1;
  639. #endif
  640.   Primitive_GC_If_Needed (Length * REAL_SIZE);
  641.   {
  642.     REAL * Temp_Array = ((REAL *) Free);
  643.     Image_Rotate_90cclw (Array, Temp_Array, nrows, ncols);
  644.     C_Array_Copy (Temp_Array, Array, Length);
  645.   }
  646.   {
  647.     SCHEME_OBJECT argument = (ARG_REF (1));
  648.     SET_PAIR_CAR (argument, (LONG_TO_UNSIGNED_FIXNUM (ncols)));
  649.     SET_PAIR_CAR ((PAIR_CDR (argument)), (LONG_TO_UNSIGNED_FIXNUM (nrows)));
  650.     PRIMITIVE_RETURN (argument);
  651.   }
  652. }
  653.  
  654. DEFINE_PRIMITIVE ("IMAGE-MIRROR!", Prim_image_mirror, 1, 1, 0)
  655. {
  656.   long nrows;
  657.   long ncols;
  658.   REAL * Array;
  659.   PRIMITIVE_HEADER (1);
  660.   {
  661.     SCHEME_OBJECT Parray;
  662.     arg_image (1, (&nrows), (&ncols), (&Parray));
  663.     Array = (ARRAY_CONTENTS (Parray));
  664.   }
  665.   C_Mirror_Image (Array, nrows, ncols);    /* side-effecting... */
  666.   PRIMITIVE_RETURN (ARG_REF (1));
  667. }
  668.  
  669.  
  670. /* C routines   referred to above  */
  671.  
  672. /* IMAGE_FAST_TRANSPOSE
  673.    A(i,j) <-> A(j,i) .
  674.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  675.    convention:= fix row & go by columns .
  676.    UNWRAP is a bijection from the compact plane to the compact interval. */
  677.  
  678. Image_Fast_Transpose (Array, nrows)       /* for square images */
  679.      REAL *Array;
  680.      long nrows;
  681. { long i, j;
  682.   long from, to;
  683.   REAL temp;
  684.   for (i=0;i<nrows;i++) {
  685.     for (j=i;j<nrows;j++) {
  686.       from = i*nrows + j;
  687.       to   = j*nrows + i;    /* (columns transposed-image) = ncols */
  688.       temp        = Array[from];
  689.       Array[from] = Array[to];
  690.       Array[to]   = temp;
  691.     }}
  692. }
  693.  
  694. /* IMAGE_TRANSPOSE
  695.    A(i,j) -> B(j,i) .
  696.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  697.    convention:= fix row & go by columns .
  698.    UNWRAP is a bijection from the compact plane to the compact interval. */
  699.  
  700. Image_Transpose (Array, New_Array, nrows, ncols)
  701.      REAL *Array, *New_Array;
  702.      long nrows, ncols;
  703. { long i, j;
  704.   for (i=0;i<nrows;i++) {
  705.     for (j=0;j<ncols;j++) {
  706.       /* (columns transposed-image) = nrows */
  707.       New_Array[j*nrows + i] = Array[i*ncols + j];
  708.     }}
  709. }
  710.  
  711. /* IMAGE_ROTATE_90CLW
  712.    A(i,j) <-> A(j, (nrows-1)-i) .
  713.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  714.    convention:= fix row & go by columns
  715.    UNWRAP is a bijection from the compact plane to the compact interval. */
  716.  
  717. Image_Rotate_90clw (Array, Rotated_Array, nrows, ncols)
  718.      REAL *Array, *Rotated_Array;
  719.      long nrows, ncols;
  720. { long i, j;
  721.  
  722.   for (i=0;i<nrows;i++) {
  723.     for (j=0;j<ncols;j++) {
  724.       /* (columns rotated_image) =nrows */
  725.       Rotated_Array[(j*nrows) + ((nrows-1)-i)] = Array[i*ncols+j];
  726.     }}
  727. }
  728.  
  729. /* ROTATION 90degrees COUNTER-CLOCK-WISE:
  730.    A(i,j) <-> A((nrows-1)-j, i) . (minus 1 because we start from 0).
  731.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  732.    because of convention:= fix row & go by columns
  733.    UNWRAP is a bijection from the compact plane to the compact interval. */
  734.  
  735. Image_Rotate_90cclw (Array, Rotated_Array, nrows, ncols)
  736.      REAL *Array, *Rotated_Array;
  737.      long nrows, ncols;
  738. { long i, j;
  739.   fast long from_index, to_index;
  740.   long Length=nrows*ncols;
  741.   for (i=0;i<nrows;i++) {
  742.     for (j=0;j<ncols;j++) {
  743.       from_index = i*ncols +j;
  744.       /* (columns rotated-image) = nrows */
  745.       to_index   = ((ncols-1)-j)*nrows + i;
  746.       Rotated_Array[to_index] = Array[from_index];
  747.     }}
  748. }
  749.  
  750. /* IMAGE_MIRROR:
  751.    A(i,j) <-> A(i, (ncols-1)-j)  [ The -1 is there because we count from 0] .
  752.    A(i,j) -------> Array[i*ncols + j]    fix row, read column convention. */
  753.  
  754. C_Mirror_Image (Array, nrows, ncols)
  755.      REAL *Array;
  756.      long nrows, ncols;
  757. { long i, j;
  758.   long ncols2=ncols/2, Length=nrows*ncols;
  759.   REAL temp;
  760.   long from, to;
  761.  
  762.   for (i=0; i<Length; i += ncols) {
  763.     for (j=0; j<ncols2; j++) {    /* DO NOT UNDO the reflections */
  764.       from = i + j;                       /* i is really i*nrows */
  765.       to   = i + (ncols-1)-j;
  766.       temp        = Array[from];
  767.       Array[from] = Array[to];
  768.       Array[to]   = temp;
  769.     }}
  770. }
  771.  
  772. /* IMAGE_ROTATE_90CLW_MIRROR:
  773.    A(i,j) <-> A(j, i)
  774.    this should be identical to image_transpose (see above).
  775.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  776.    because of convention:= fix row & go by columns
  777.    UNWRAP is a bijection from the compact plane to the compact interval. */
  778.  
  779. C_Rotate_90clw_Mirror_Image (Array, Rotated_Array, nrows, ncols)
  780.      REAL *Array, *Rotated_Array;
  781.      long nrows, ncols;
  782. { long i, j;
  783.   long from, to, Length=nrows*ncols;
  784.  
  785.   for (i=0;i<nrows;i++) {
  786.     for (j=0;j<ncols;j++) {
  787.       from = i*ncols +j;
  788.       /* the columns of the rotated image are nrows! */
  789.       to   = j*nrows +i;
  790.       Rotated_Array[to] = Array[from];
  791.     }}
  792. }
  793.  
  794. DEFINE_PRIMITIVE ("SQUARE-IMAGE-TIME-REVERSE!", Prim_square_image_time_reverse, 2,2, 0)
  795. { long i, rows;
  796.   REAL *a;
  797.   void square_image_time_reverse();
  798.   PRIMITIVE_HEADER (2);
  799.   CHECK_ARG (1, ARRAY_P);
  800.   CHECK_ARG (2, FIXNUM_P);
  801.   a = ARRAY_CONTENTS(ARG_REF(1));
  802.   rows = arg_nonnegative_integer(2);
  803.   if ((rows*rows) != ARRAY_LENGTH(ARG_REF(1)))     error_bad_range_arg(1);
  804.   square_image_time_reverse(a,rows);
  805.  
  806.   PRIMITIVE_RETURN (UNSPECIFIC);
  807. }
  808.  
  809. /* Square Image Time reverse
  810.    is combination of one-dimensional time-reverse
  811.    row-wise and column-wise.
  812.    It can be done slightly more efficiently than below. */
  813.  
  814. void
  815. square_image_time_reverse (x,rows)
  816.      REAL *x;
  817.      long rows;
  818. { long i,cols;
  819.   REAL *xrow, *yrow;
  820.   void C_Array_Time_Reverse();
  821.   cols = rows;            /* square image */
  822.  
  823.   xrow = x;
  824.   for (i=0; i<rows; i++)    /* row-wise */
  825.   { C_Array_Time_Reverse(xrow,cols);
  826.     xrow = xrow + cols; }
  827.  
  828.   Image_Fast_Transpose(x, rows);
  829.  
  830.   xrow = x;
  831.   for (i=0; i<rows; i++)    /* column-wise */
  832.   { C_Array_Time_Reverse(xrow,cols);
  833.     xrow = xrow + cols; }
  834.  
  835.   Image_Fast_Transpose(x, rows);
  836. }
  837.  
  838. /* cs-images */
  839.  
  840. /* operation-1
  841.    groups together procedures     that operate on 1 cs-image-array
  842.    (side-effecting the image) */
  843.  
  844. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-1!", Prim_cs_image_operation_1, 3,3, 0)
  845. { long rows, opcode;
  846.   REAL *a;
  847.   void cs_image_magnitude(), cs_image_real_part(), cs_image_imag_part();
  848.   PRIMITIVE_HEADER (3);
  849.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  850.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  851.   CHECK_ARG (3, ARRAY_P);    /* input and output image array */
  852.  
  853.   a = ARRAY_CONTENTS(ARG_REF(3));
  854.   rows = arg_nonnegative_integer(2); /*          square images only */
  855.   if ((rows*rows) != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(1);
  856.   opcode = arg_nonnegative_integer(1);
  857.  
  858.   if (opcode==1)
  859.     cs_image_magnitude(a,rows);
  860.   else if (opcode==2)
  861.     cs_image_real_part(a,rows);
  862.   else if (opcode==3)
  863.     cs_image_imag_part(a,rows);
  864.   else
  865.     error_bad_range_arg(3);    /* illegal opcode */
  866.  
  867.   PRIMITIVE_RETURN (UNSPECIFIC);
  868. }
  869.  
  870.  
  871. void
  872. cs_array_real_part (a,n)
  873.      REAL *a;
  874.      long n;
  875. { long i,n2;            /* works both for even and odd length */
  876.   n2 = n/2;
  877.   for (i=n2+1;i<n;i++) a[i] = a[n-i]; /* copy real values into place */
  878.   /*                                     even signal */
  879. }
  880.  
  881. void
  882. cs_array_imag_part (a,n)
  883.      REAL *a;
  884.      long n;
  885. { long i,n2;
  886.   n2 = n/2;            /* integer division truncates down */
  887.   for (i=n2+1; i<n; i++)    /* works both for even and odd length */
  888.   { a[n-i] = a[i];        /* copy imaginary values into place */
  889.     a[i]   = (-a[i]); }        /* odd signal */
  890.   a[0]     = 0.0;
  891.   if (2*n2 == n)        /* even length, n2 is real only */
  892.     a[n2]    = 0.0;
  893. }
  894.  
  895. /* From now on (below), assume that cs-images (rows=cols)
  896.    have always EVEN LENGTH which is true when they come from FFTs */
  897.  
  898. /*  In the following 3   time-reverse the bottom half rows
  899.     is done to match  the frequencies of complex-images
  900.     coming from cft2d.
  901.     Also transpose is needed to match frequencies identically
  902.  
  903.     #|
  904.     ;; Scrabling of frequencies in  cs-images
  905.  
  906.     ;; start from real image  4x4
  907.  
  908.     ;; rft2d    is a cs-image
  909.     (3.5 .375 -2.75 1.875    -.25 0. 0. -.25    -.25 -.125 0. .125
  910.     .25 .25 0. 0.)
  911.  
  912.     ;; cft2d   transposed
  913.     ;; real
  914.     3.5 .375 -2.75 .375
  915.     -.25  0.  0.  -.25  ; same as cs-image
  916.     -.25 -.125 0. -.125
  917.     -.25 -.25  0.   0.  ; row3 = copy 1 + time-reverse
  918.     ;; imag
  919.     0. 1.875 0. -1.875
  920.     .25 .25 0. 0.       ; same as cs-image
  921.     0. .125 0. -.125
  922.     -.25 0. 0. -.25     ; row 3 = copy 1 + negate + time-reverse
  923.     |#
  924.  
  925.     */
  926.  
  927. void
  928. cs_image_magnitude (x,rows)
  929.      REAL *x;
  930.      long rows;
  931. { long i,j, cols, n,n2, nj; /*     result = real ordinary image */
  932.   REAL *xrow, *yrow;
  933.   cols = rows;            /* input cs-image   is square */
  934.   n = rows;
  935.   n2 = n/2;
  936.  
  937.   xrow = x;
  938.   cs_array_magnitude(xrow, n);  /* row 0 is cs-array */
  939.   xrow = x + n2*cols;
  940.   cs_array_magnitude(xrow, n);  /* row n2 is cs-array */
  941.  
  942.   xrow = x + cols;        /* real part */
  943.   yrow = x + (rows-1)*cols;    /* imag part */
  944.   for (i=1; i<n2; i++) {
  945.     xrow[ 0] = (REAL) sqrt((double) xrow[ 0]*xrow[ 0] + yrow[ 0]*yrow[ 0]);
  946.     xrow[n2] = (REAL) sqrt((double) xrow[n2]*xrow[n2] + yrow[n2]*yrow[n2]);
  947.     yrow[ 0] = xrow[ 0];
  948.     yrow[n2] = xrow[n2];
  949.     for (j=1; j<n2; j++) {
  950.       nj = n-j;
  951.       xrow[ j] = (REAL) sqrt((double) xrow[ j]*xrow[ j] + yrow[ j]*yrow[ j]);
  952.       xrow[nj] = (REAL) sqrt((double) xrow[nj]*xrow[nj] + yrow[nj]*yrow[nj]);
  953.       yrow[j]  = xrow[nj];
  954.       yrow[nj] = xrow[ j];      /* Bottom rows: copy (even) and time-reverse */
  955.     }
  956.     xrow = xrow + cols;
  957.     yrow = yrow - cols; }
  958.   Image_Fast_Transpose(x, n);
  959. }
  960.  
  961. void
  962. cs_image_real_part (x,rows)
  963.      REAL *x;
  964.      long rows;
  965. { long i,j,cols, n,n2;
  966.   REAL *xrow, *yrow;
  967.   void cs_array_real_part();
  968.   cols = rows;            /* square image */
  969.   n = rows;
  970.   n2 = n/2;
  971.  
  972.   xrow = x;
  973.   cs_array_real_part(xrow, n);  /* row 0 is cs-array */
  974.   xrow = x + n2*cols;
  975.   cs_array_real_part(xrow, n);  /* row n2 is cs-array */
  976.  
  977.   xrow = x + cols;        /* real part */
  978.   yrow = x + (rows-1)*cols;    /* imag part */
  979.   for (i=1; i<n2; i++) {
  980.     /* copy real part into imaginary's place (even) */
  981.     yrow[0]  = xrow[0];
  982.     for (j=1; j<n; j++)
  983.       yrow[j] = xrow[n-j];    /* Bottom rows:  copy and time-reverse */
  984.     xrow = xrow + cols;
  985.     yrow = yrow - cols; }
  986.   Image_Fast_Transpose(x, n);
  987. }
  988.  
  989. void
  990. cs_image_imag_part (x,rows)
  991.      REAL *x;
  992.      long rows;
  993. { long i,j,cols, n,n2, nj;
  994.   REAL *xrow, *yrow;
  995.   void cs_array_imag_part();
  996.   cols = rows;            /* square image */
  997.   n = rows;
  998.   n2 = n/2;
  999.  
  1000.   xrow = x;
  1001.   cs_array_imag_part(xrow, n);  /* row 0 is cs-array */
  1002.   xrow = x + n2*cols;
  1003.   cs_array_imag_part(xrow, n);  /* row n2 is cs-array */
  1004.  
  1005.   xrow = x + cols;        /* real part */
  1006.   yrow = x + (rows-1)*cols;    /* imag part */
  1007.   for (i=1; i<n2; i++) {
  1008.     xrow[0]  = yrow[0];        /* copy the imaginary part into real's place */
  1009.     xrow[n2] = yrow[n2];
  1010.     yrow[0]  = (-yrow[0]);      /* negate (odd) */
  1011.     yrow[n2] = (-yrow[n2]);
  1012.     for (j=1;j<n2; j++) {
  1013.       nj = n-j;
  1014.       xrow[j]  = yrow[j];     /* copy the imaginary part into real's place */
  1015.       xrow[nj] = yrow[nj];
  1016.       /* Bottom rows: negate (odd) and time-reverse */
  1017.       yrow[j]  = (-xrow[nj]);
  1018.       yrow[nj] = (-xrow[j]); }
  1019.     xrow = xrow + cols;
  1020.     yrow = yrow - cols; }
  1021.   Image_Fast_Transpose(x, n);
  1022. }
  1023.  
  1024. /* cs-image-operation-2
  1025.    groups together procedures     that use 2 cs-image-arrays
  1026.    (usually side-effecting the 2nd image, but not necessarily) */
  1027.  
  1028. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2!", Prim_cs_image_operation_2, 4, 4, 0)
  1029. { long rows, nn, opcode;
  1030.   REAL *x,*y;
  1031.   void cs_image_multiply_into_second_one();
  1032.   PRIMITIVE_HEADER (4);
  1033.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1034.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1035.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1036.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1037.  
  1038.   x = ARRAY_CONTENTS(ARG_REF(3));
  1039.   y = ARRAY_CONTENTS(ARG_REF(4));
  1040.   rows = arg_nonnegative_integer(2); /*          square images only */
  1041.   nn = rows*rows;
  1042.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1043.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1044.  
  1045.   opcode = arg_nonnegative_integer(1);
  1046.  
  1047.   if (opcode==1)
  1048.     cs_image_multiply_into_second_one(x,y,rows); /* result in y */
  1049.   else if (opcode==2)
  1050.     error_bad_range_arg(1);    /* illegal opcode */
  1051.   else
  1052.     error_bad_range_arg(1);    /* illegal opcode */
  1053.  
  1054.   PRIMITIVE_RETURN (UNSPECIFIC);
  1055. }
  1056.  
  1057. void
  1058. cs_image_multiply_into_second_one (x,y, rows)
  1059.      REAL *x,*y;
  1060.      long rows;
  1061. { long i,j,cols, n,n2;
  1062.   REAL *xrow,*yrow,  *xrow_r, *xrow_i, *yrow_r, *yrow_i, temp;
  1063.   cols = rows;            /* square image */
  1064.   n = rows;
  1065.   n2 = n/2;
  1066.  
  1067.   xrow= x; yrow= y;
  1068.   cs_array_multiply_into_second_one(xrow,yrow, n,n2); /*         row 0 */
  1069.  
  1070.   xrow= x+n2*cols; yrow= y+n2*cols;
  1071.   cs_array_multiply_into_second_one(xrow,yrow, n,n2); /*         row n2 */
  1072.  
  1073.   xrow_r= x+cols;           yrow_r= y+cols;
  1074.   xrow_i= x+(n-1)*cols;     yrow_i= y+(n-1)*cols;
  1075.   for (i=1; i<n2; i++) {
  1076.     for (j=0; j<n; j++) {
  1077.       temp      = xrow_r[j]*yrow_r[j]  -  xrow_i[j]*yrow_i[j]; /* real part */
  1078.       yrow_i[j] = xrow_r[j]*yrow_i[j]  +  xrow_i[j]*yrow_r[j]; /* imag part */
  1079.       yrow_r[j] = temp; }
  1080.     xrow_r= xrow_r+cols;   yrow_r= yrow_r+cols;
  1081.     xrow_i= xrow_i-cols;   yrow_i= yrow_i-cols;
  1082.   }
  1083. }
  1084.  
  1085. /* cs-image-operation-2x!     is just like     cs-image-operation-2!
  1086.   but takes an additional flonum argument. */
  1087.  
  1088. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2x!", Prim_cs_image_operation_2x, 5, 5, 0)
  1089. { long rows, nn, opcode;
  1090.   REAL *x,*y, flonum_arg;
  1091.   void cs_image_divide_into_z();
  1092.   PRIMITIVE_HEADER (5);
  1093.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1094.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1095.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1096.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1097.   flonum_arg = (arg_real (5));
  1098.  
  1099.   x = ARRAY_CONTENTS(ARG_REF(3));
  1100.   y = ARRAY_CONTENTS(ARG_REF(4));
  1101.   rows = arg_nonnegative_integer(2); /*          square images only */
  1102.   nn = rows*rows;
  1103.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1104.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1105.  
  1106.   opcode = arg_nonnegative_integer(1);
  1107.  
  1108.   if (opcode==1)
  1109.     cs_image_divide_into_z( x,y, x, rows, flonum_arg); /* result in x */
  1110.   else if (opcode==2)
  1111.     cs_image_divide_into_z( x,y, y, rows, flonum_arg); /* result in y */
  1112.   else
  1113.     error_bad_range_arg(1);    /* illegal opcode */
  1114.  
  1115.   PRIMITIVE_RETURN (UNSPECIFIC);
  1116. }
  1117.  
  1118. /* The convention for inf values in division 1/0 is just like in arrays */
  1119.  
  1120. void
  1121. cs_image_divide_into_z (x,y, z, rows, inf)
  1122.      REAL *x,*y,*z, inf;    /* z can be either x or y */
  1123.      long rows;
  1124. { long i,j,cols, n,n2;
  1125.   REAL temp, radius;
  1126.   REAL  *ar_,*ai_, *br_,*bi_, *zr_,*zi_; /* Letters a,b  correspond to  x,y */
  1127.   REAL *xrow,*yrow,*zrow;
  1128.   cols = rows;            /* square image */
  1129.   n = rows;
  1130.   n2 = n/2;
  1131.  
  1132.   xrow= x; yrow= y; zrow= z;
  1133.   cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /*         row 0 */
  1134.  
  1135.   xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
  1136.   cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /*         row n2 */
  1137.  
  1138.   ar_= x+cols;           br_= y+cols;            zr_= z+cols;
  1139.   ai_= x+(n-1)*cols;     bi_= y+(n-1)*cols;      zi_= z+(n-1)*cols;
  1140.   for (i=1; i<n2; i++) {
  1141.     for (j=0; j<n; j++) {
  1142.       /* b^2 denominator = real^2 + imag^2 */
  1143.       radius = br_[j]*br_[j]  + bi_[j]*bi_[j];
  1144.  
  1145.       if (radius == 0.0) {
  1146.     if (ar_[j] == 0.0)  zr_[j]  = 1.0;
  1147.     else                zr_[j]  = ar_[j] * inf;
  1148.     if (ai_[j] == 0.0)  zi_[j]  = 1.0;
  1149.     else                zi_[j]  = ai_[j] * inf; }
  1150.       else {
  1151.     temp    =  ar_[j]*br_[j]   +  ai_[j]*bi_[j];
  1152.     zi_[j]  = (ai_[j]*br_[j]   -  ar_[j]*bi_[j]) / radius; /* imag part */
  1153.     zr_[j]  = temp                               / radius; /* real part */
  1154.       }}
  1155.     ar_= ar_+cols;   br_= br_+cols;    zr_= zr_+cols;
  1156.     ai_= ai_-cols;   bi_= bi_-cols;    zi_= zi_-cols;
  1157.   }
  1158. }
  1159.  
  1160. /* operation-3
  1161.    groups together procedures     that use 3 cs-image-arrays
  1162.    (usually side-effecting the 3rd image, but not necessarily) */
  1163.  
  1164. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-3!", Prim_cs_image_operation_3, 5, 5, 0)
  1165. { long rows, nn, opcode;
  1166.   REAL *x,*y,*z;
  1167.   void tr_complex_image_to_cs_image();
  1168.   PRIMITIVE_HEADER (5);
  1169.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1170.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1171.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1172.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1173.   CHECK_ARG (5, ARRAY_P);    /* image array 3 */
  1174.  
  1175.   x = ARRAY_CONTENTS(ARG_REF(3));
  1176.   y = ARRAY_CONTENTS(ARG_REF(4));
  1177.   z = ARRAY_CONTENTS(ARG_REF(5));
  1178.   rows = arg_nonnegative_integer(2); /*          square images only */
  1179.   nn = rows*rows;
  1180.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1181.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1182.   if (nn != ARRAY_LENGTH(ARG_REF(5)))   error_bad_range_arg(5);
  1183.  
  1184.   opcode = arg_nonnegative_integer(1);
  1185.  
  1186.   if (opcode==1)
  1187.     tr_complex_image_to_cs_image(x,y, z,rows); /* result in z */
  1188.   else if (opcode==2)
  1189.     error_bad_range_arg(1);    /* illegal opcode */
  1190.   else
  1191.     error_bad_range_arg(1);    /* illegal opcode */
  1192.  
  1193.   PRIMITIVE_RETURN (UNSPECIFIC);
  1194. }
  1195.  
  1196. /* x and y must be ALREADY TRANSPOSED real and imaginary parts */
  1197. void
  1198. tr_complex_image_to_cs_image (x,y, z,rows)
  1199.      REAL *x,*y,*z;
  1200.      long rows;
  1201. { long i,j,cols, n,n2, n2_1_n;
  1202.   REAL *xrow, *yrow, *zrow;
  1203.   cols = rows;            /* square image */
  1204.   n = rows;
  1205.   n2 = n/2;
  1206.  
  1207.   xrow= x; yrow= y; zrow= z;
  1208.   for (j=0; j<=n2; j++)
  1209.     /* real part of row 0 (cs-array) */
  1210.     zrow[j] = xrow[j];
  1211.   for (j=n2+1; j<n; j++)
  1212.     /* imag part of row 0 */
  1213.     zrow[j] = yrow[n-j];
  1214.   xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
  1215.   for (j=0; j<=n2; j++)
  1216.     /* real part of row n2 (cs-array) */
  1217.     zrow[j] = xrow[j];
  1218.   for (j=n2+1; j<n; j++)
  1219.     /* imag part of row n2 */
  1220.     zrow[j] = yrow[n-j];
  1221.   xrow= x+cols;   zrow= z+cols;   n2_1_n = (n2-1)*cols;
  1222.   for (j=0; j<n2_1_n; j++)
  1223.     /* real rows 1,2,..,n2-1 */
  1224.     zrow[j] = xrow[j];
  1225.   yrow= y+(n2-1)*cols;
  1226.   /* imag rows n2+1,n2+2,... */
  1227.   zrow= z+(n2+1)*cols;
  1228.   for (i=1; i<n2; i++) {
  1229.     for (j=0; j<n; j++)   zrow[j] = yrow[j];
  1230.     zrow = zrow + cols;
  1231.     yrow = yrow - cols;
  1232.   }
  1233. }
  1234.